#ifndef SPLINES_H
#define SPLINES_H

#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <vector>
#include <fstream>
#include <math.h>
using namespace std;

class Splines;
class Function;
double fB1(vector<double> v, int i, double x);
double fB2(vector<double> v, int i, double x);
double dfB2(vector<double> v, int i, double x);
double fB3(vector<double> v, int i, double x);
double dfB3(vector<double> v, int i, double x);
double ddfB3(vector<double> v, int i, double x);
Splines interpolation(vector<double> x, vector<double> y, double dl, double dr, double ddl, double ddr, string algorithm, string spline, string boundary);
vector<double> solve1(vector<double> a, vector<double> b);
vector<double> solve2(double x1, double x2, double f1, double f2, double m1, double m2);
vector<double> solve3(double x1, double x2, double f1, double f2, double M1, double M2);
void pushback(vector<double> &x, double a[], int n);
void printvector(vector<double> v);

//----------------------------------------------------------------------------------------------------------//

// 样条
class Splines
{
public:
    string spline;                  // linear or cubic
    vector<double> x;               // vector x
    vector<vector<double>> coe;     // coefficient vector
    void printspline();             // print coefficient vectors
    void fileprintspline(string s); // print n,x,coe into files
    double value(double x);         // derive the value of x
};

// 函数
class Function
{
public:
    virtual double f(double x)
    {
        return 0;
    }
    virtual double df(double x)
    {
        return 0;
    }
    virtual double ddf(double x)
    {
        return 0;
    }
};

// B1 values
double fB1(vector<double> v, int i, double x)
{
    int n = v.size();
    if (i < 0 || i > n - 1)
        return 0;
    else if (i == 0)
    {
        if (x > v[0] && x <= v[1])
            return (v[1] - x) / (v[1] - v[0]);
        else
            return 0;
    }
    else if (i == n - 1)
    {
        if (x > v[n - 2] && x <= v[n - 1])
            return (x - v[n - 2]) / (v[n - 1] - v[n - 2]);
        else
            return 0;
    }
    else
    {
        if (x > v[i - 1] && x <= v[i])
            return (x - v[i - 1]) / (v[i] - v[i - 1]);
        else if (x > v[i] && x <= v[i + 1])
            return (v[i + 1] - x) / (v[i + 1] - v[x]);
        else
            return 0;
    }
}
// B2 values and derivatives
double fB2(vector<double> v, int i, double x)
{
    int n = v.size();
    if (i < 0 || i >= n - 1)
        return 0;
    else if (i == 0)
        return fB1(v, 1, x) * (v[2] - x) / (v[2] - v[0]);
    else if (i == n - 2)
        return fB1(v, n - 2, x) * (x - v[n - 3]) / (v[n - 1] - v[n - 3]);
    else
        return fB1(v, i, x) * (x - v[i - 1]) / (v[i + 1] - v[i - 1]) + fB1(v, i + 1, x) * (v[i + 2] - x) / (v[i + 2] - v[i]);
}
double dfB2(vector<double> v, int i, double x)
{
    int n = v.size();
    if (i < 0 || i >= n - 1)
        return 0;
    else if (i == 0)
        return -2 * fB1(v, 1, x) / (v[2] - v[0]);
    else if (i == n - 2)
        return 2 * fB1(v, n - 2, x) / (v[n - 1] - v[n - 3]);
    else
        return 2 * fB1(v, i, x) / (v[i + 1] - v[i - 1]) - 2 * fB1(v, i + 1, x) / (v[i + 2] - v[i]);
}
// B3 values and derivatives
double fB3(vector<double> v, int i, double x)
{
    int n = v.size();
    if (i < 0 || i >= n - 2)
        return 0;
    else if (i == 0)
        return fB2(v, 1, x) * (v[3] - x) / (v[3] - v[0]);
    else if (i == n - 3)
        return fB2(v, n - 3, x) * (x - v[n - 4]) / (v[n - 1] - v[n - 4]);
    else
        return fB2(v, i, x) * (x - v[i - 1]) / (v[i + 2] - v[i - 1]) + fB2(v, i + 1, x) * (v[i + 3] - x) / (v[i + 3] - v[i]);
}
double dfB3(vector<double> v, int i, double x)
{
    int n = v.size();
    if (i < 0 || i >= n - 2)
        return 0;
    else if (i == 0)
        return -3 * fB2(v, 1, x) / (v[3] - v[0]);
    else if (i == n - 3)
        return 3 * fB2(v, n - 3, x) / (v[n - 1] - v[n - 4]);
    else
        return 3 * fB2(v, i, x) / (v[i + 2] - v[i - 1]) - 3 * fB2(v, i + 1, x) / (v[i + 3] - v[i]);
}
double ddfB3(vector<double> v, int i, double x)
{
    int n = v.size();
    if (i < 0 || i >= n - 2)
        return 0;
    else if (i == 0)
        return -3 * dfB2(v, 1, x) / (v[3] - v[0]);
    else if (i == n - 3)
        return 3 * dfB2(v, n - 3, x) / (v[n - 1] - v[n - 4]);
    else
        return 3 * dfB2(v, i, x) / (v[i + 2] - v[i - 1]) - 3 * dfB2(v, i + 1, x) / (v[i + 3] - v[i]);
}

// interpolation
// algorithm: ppForm or Bsplines
// spline: linear or cubic
// boundary: ccs(complete cubic spline) or csssd(cubic spline with specified second derivatives)
//  or ncs(natural cubic spline)
Splines interpolation(vector<double> x, vector<double> y, double dl, double dr, double ddl, double ddr, string algorithm, string spline, string boundary)
{
    Splines s;
    s.x = x;
    s.spline = spline;
    int n = x.size();
    vector<double> v; // 系数
    if (algorithm == "ppForm")
    {
        if (spline == "linear")
        {
            double c1, c0; // 一次项和常数项系数
            for (int i = 0; i <= n - 2; i++)
            {
                c1 = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
                c0 = y[i] - c1 * x[i];
                v.push_back(c1);
                v.push_back(c0);
                s.coe.push_back(v);
                v.clear();
            }
        }
        else if (spline == "cubic")
        {
            vector<double> p(n - 1), q(n - 1), a, b, df(n - 1), ddf(n - 2);
            for (int i = 1; i <= n - 2; i++)
            {
                p[i] = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]); // miu
                q[i] = (x[i + 1] - x[i]) / (x[i + 1] - x[i - 1]); // lamda
            }
            for (int i = 0; i <= n - 2; i++)
                df[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]); // f([x_i,x_i+1])
            for (int i = 0; i <= n - 3; i++)
                ddf[i] = (df[i + 1] - df[i]) / (x[i + 2] - x[i]); // f([x_i,x_i+1,x_i+2])
            if (boundary == "ccs")
            {
                for (int i = 1; i <= n - 3; i++)
                    a.push_back(p[i]);
                for (int i = 1; i <= n - 2; i++)
                    a.push_back(2.0);
                for (int i = 2; i <= n - 2; i++)
                    a.push_back(q[i]);
                b.push_back(3 * p[1] * df[1] + 3 * q[1] * df[0] - q[1] * dl);
                for (int i = 2; i <= n - 3; i++)
                    b.push_back(3 * p[i] * df[i] + 3 * q[i] * df[i - 1]);
                b.push_back(3 * p[n - 2] * df[n - 2] + 3 * q[n - 2] * df[n - 3] - p[n - 2] * dr);
                vector<double> c, d;
                c = solve1(a, b);
                d.push_back(dl);
                for (int i = 1; i <= n - 2; i++)
                    d.push_back(c[i - 1]);
                d.push_back(dr);
                for (int i = 0; i <= n - 2; i++)
                {
                    v = solve2(x[i], x[i + 1], y[i], y[i + 1], d[i], d[i + 1]);
                    s.coe.push_back(v);
                    v.clear();
                }
            }
            else if (boundary == "csssd" || boundary == "ncs")
            {
                for (int i = 1; i <= n - 3; i++)
                    a.push_back(q[i]);
                for (int i = 1; i <= n - 2; i++)
                    a.push_back(2.0);
                for (int i = 2; i <= n - 2; i++)
                    a.push_back(p[i]);
                if (boundary == "csssd")
                    b.push_back(6 * ddf[0] - p[1] * ddl);
                else
                    b.push_back(6 * ddf[0]);
                for (int i = 1; i <= n - 4; i++)
                    b.push_back(6 * ddf[i]);
                if (boundary == "csssd")
                    b.push_back(6 * ddf[n - 3] - q[n - 2] * ddr);
                else
                    b.push_back(6 * ddf[n - 3]);
                vector<double> c, d;
                c = solve1(a, b);
                if (boundary == "csssd")
                    d.push_back(ddl);
                else
                    d.push_back(0);
                for (int i = 1; i <= n - 2; i++)
                    d.push_back(c[i - 1]);
                if (boundary == "csssd")
                    d.push_back(ddr);
                else
                    d.push_back(0);
                for (int i = 0; i <= n - 2; i++)
                {
                    v = solve3(x[i], x[i + 1], y[i], y[i + 1], d[i], d[i + 1]);
                    s.coe.push_back(v);
                    v.clear();
                }
            }
        }
    }
    else if (algorithm == "Bsplines")
    {
        if (spline == "linear")
        {
            // a_i=f_i
            double c1, c0; // 一次项和常数项系数
            for (int i = 0; i <= n - 2; i++)
            {
                c1 = y[i + 1] / (x[i + 1] - x[i]) - y[i] / (x[i + 1] - x[i]);
                c0 = y[i] * x[i + 1] / (x[i + 1] - x[i]) - y[i + 1] * x[i] / (x[i + 1] - x[i]);
                v.push_back(c1);
                v.push_back(c0);
                s.coe.push_back(v);
                v.clear();
            }
        }
        else if (spline == "cubic")
        {
            vector<double> a, b, c, d;
            a.push_back(x[0] - 2.1), a.push_back(x[0] - 1);
            for (int i = 0; i < n; i++)
                a.push_back(x[i]);
            a.push_back(x[n - 1] + 1), a.push_back(x[n - 1] + 2.1);
            b.push_back(0);
            for (int i = 0; i < n; i++)
                b.push_back(y[i]);
            b.push_back(0);
            if (boundary == "ccs")
            {
                c.push_back(dfB3(a, 1, x[0]) - dfB3(a, 2, x[0]) * fB3(a, 1, x[0]) / fB3(a, 2, x[0]));
                for (int i = 0; i < n; i++)
                    c.push_back(fB3(a, i + 2, x[i]));
                c.push_back(dfB3(a, 0, x[0]) - dfB3(a, 2, x[0]) * fB3(a, 0, x[0]) / fB3(a, 2, x[0]));
                for (int i = 0; i < n; i++)
                    c.push_back(fB3(a, i + 1, x[i]));
                c.push_back(dfB3(a, n + 1, x[n - 1]) - dfB3(a, n - 1, x[n - 1]) * fB3(a, n + 1, x[n - 1]) / fB3(a, n - 1, x[n - 1]));
                for (int i = 0; i < n; i++)
                    c.push_back(fB3(a, i, x[i]));
                c.push_back(dfB3(a, n, x[n - 1]) - dfB3(a, n - 1, x[n - 1]) * fB3(a, n, x[n - 1]) / fB3(a, n - 1, x[n - 1]));
                b[0] = dl - y[0] * dfB3(a, 2, x[0]) / fB3(a, 2, x[0]);
                b[n + 1] = dr - y[n - 1] * dfB3(a, n - 1, x[n - 1]) / fB3(a, n - 1, x[n - 1]);
                d = solve1(c, b);
            }
            else if (boundary == "csssd" || boundary == "ncs")
            {
                c.push_back(ddfB3(a, 1, x[0]) - ddfB3(a, 2, x[0]) * fB3(a, 1, x[0]) / fB3(a, 2, x[0]));
                for (int i = 0; i < n; i++)
                    c.push_back(fB3(a, i + 2, x[i]));
                c.push_back(ddfB3(a, 0, x[0]) - ddfB3(a, 2, x[0]) * fB3(a, 0, x[0]) / fB3(a, 2, x[0]));
                for (int i = 0; i < n; i++)
                    c.push_back(fB3(a, i + 1, x[i]));
                c.push_back(ddfB3(a, n + 1, x[n - 1]) - ddfB3(a, n - 1, x[n - 1]) * fB3(a, n + 1, x[n - 1]) / fB3(a, n - 1, x[n - 1]));
                for (int i = 0; i < n; i++)
                    c.push_back(fB3(a, i, x[i]));
                c.push_back(ddfB3(a, n, x[n - 1]) - ddfB3(a, n - 1, x[n - 1]) * fB3(a, n, x[n - 1]) / fB3(a, n - 1, x[n - 1]));
                if (boundary == "csssd")
                {
                    b[0] = ddl - y[0] * ddfB3(a, 2, x[0]) / fB3(a, 2, x[0]);
                    b[n + 1] = ddr - y[n - 1] * ddfB3(a, n - 1, x[n - 1]) / fB3(a, n - 1, x[n - 1]);
                }
                else
                {
                    b[0] = -y[0] * ddfB3(a, 2, x[0]) / fB3(a, 2, x[0]);
                    b[n + 1] = -y[n - 1] * ddfB3(a, n - 1, x[n - 1]) / fB3(a, n - 1, x[n - 1]);
                }
                d = solve1(c, b);
            }
            for (int i = 2; i <= n; i++)
            {
                vector<double> v1(4), v2(4), v3(4), v4(4);
                double k;
                k = (a[i + 1] - a[i - 2]) * (a[i + 1] - a[i - 1]) * (a[i + 1] - a[i]);
                v1[0] = -1 / k, v1[1] = 3 * a[i + 1] / k, v1[2] = -3 * a[i + 1] * a[i + 1] / k, v1[3] = a[i + 1] * a[i + 1] * a[i + 1] / k;
                k = (a[i + 1] - a[i - 2]) * (a[i + 1] - a[i - 1]) * (a[i + 1] - a[i]);
                v2[0] += 1 / k, v2[1] -= (a[i - 2] + 2 * a[i + 1]) / k, v2[2] += a[i + 1] * (a[i + 1] + 2 * a[i - 2]) / k,
                    v2[3] -= a[i - 2] * a[i + 1] * a[i + 1] / k;
                k = (a[i + 2] - a[i - 1]) * (a[i + 1] - a[i - 1]) * (a[i + 1] - a[i]);
                v2[0] += 1 / k, v2[1] -= (a[i + 2] + a[i + 1] + a[i - 1]) / k, v2[2] += (a[i + 2] * a[i + 1] + a[i + 2] * a[i - 1] + a[i + 1] * a[i - 1]) / k,
                    v2[3] -= a[i + 2] * a[i + 1] * a[i - 1] / k;
                k = (a[i + 2] - a[i]) * (a[i + 1] - a[i]) * (a[i + 2] - a[i - 1]);
                v2[0] += 1 / k, v2[1] -= (a[i] + 2 * a[i + 2]) / k, v2[2] += a[i + 2] * (a[i + 2] + 2 * a[i]) / k,
                    v2[3] -= a[i] * a[i + 2] * a[i + 2] / k;
                k = (a[i + 2] - a[i - 1]) * (a[i + 1] - a[i - 1]) * (a[i + 1] - a[i]);
                v3[0] -= 1 / k, v3[1] += (a[i + 1] + 2 * a[i - 1]) / k, v3[2] -= a[i - 1] * (a[i - 1] + 2 * a[i + 1]) / k,
                    v3[3] += a[i + 1] * a[i - 1] * a[i - 1] / k;
                k = (a[i + 2] - a[i - 1]) * (a[i + 2] - a[i]) * (a[i + 1] - a[i]);
                v3[0] -= 1 / k, v3[1] += (a[i + 2] + a[i] + a[i - 1]) / k, v3[2] -= (a[i + 2] * a[i] + a[i + 2] * a[i - 1] + a[i] * a[i - 1]) / k,
                    v3[3] += a[i + 2] * a[i] * a[i - 1] / k;
                k = (a[i + 3] - a[i]) * (a[i + 2] - a[i]) * (a[i + 1] - a[i]);
                v3[0] -= 1 / k, v3[1] += (a[i + 3] + 2 * a[i]) / k, v3[2] -= a[i] * (a[i] + 2 * a[i + 3]) / k,
                    v3[3] += a[i + 3] * a[i] * a[i] / k;
                k = (a[i + 3] - a[i]) * (a[i + 2] - a[i]) * (a[i + 1] - a[i]);
                v4[0] = 1 / k, v4[1] = -3 * a[i] / k, v4[2] = 3 * a[i] * a[i] / k, v4[3] = -a[i] * a[i] * a[i] / k;
                v.push_back(d[i - 2] * v1[0] + d[i - 1] * v2[0] + d[i] * v3[0] + d[i + 1] * v4[0]);
                v.push_back(d[i - 2] * v1[1] + d[i - 1] * v2[1] + d[i] * v3[1] + d[i + 1] * v4[1]);
                v.push_back(d[i - 2] * v1[2] + d[i - 1] * v2[2] + d[i] * v3[2] + d[i + 1] * v4[2]);
                v.push_back(d[i - 2] * v1[3] + d[i - 1] * v2[3] + d[i] * v3[3] + d[i + 1] * v4[3]);
                s.coe.push_back(v);
                v.clear();
            }
        }
    }
    return s;
}

// AX=b求解，其中A为三对角矩阵
vector<double> solve1(vector<double> a, vector<double> b)
{
    int n = b.size();
    vector<double> v(n), w(n), x(n), y(n), z(n), p(n), q(n), r(n);
    for (int i = 1; i <= n - 1; i++)
        x[i] = a[i + 2 * n - 2];
    for (int i = 0; i <= n - 1; i++)
        y[i] = a[i + n - 1];
    for (int i = 0; i <= n - 2; i++)
        z[i] = a[i];
    p[0] = y[0], q[0] = z[0] / y[0];
    for (int i = 1; i <= n - 1; i++)
    {
        r[i] = x[i];
        p[i] = y[i] - x[i] * q[i - 1];
        q[i] = z[i] / p[i];
    }
    w[0] = b[0] / y[0];
    for (int i = 1; i <= n - 1; i++)
        w[i] = (b[i] - x[i] * w[i - 1]) / (y[i] - x[i] * q[i - 1]);
    v[n - 1] = w[n - 1];
    for (int i = n - 2; i >= 0; i--)
        v[i] = w[i] - q[i] * v[i + 1];
    return v;
}

// ppForm中可求m_i的情况求系数
vector<double> solve2(double x1, double x2, double f1, double f2, double m1, double m2)
{
    vector<double> v;
    double a, b, c, d;
    double k = (f1 - f2) / (x1 - x2);
    a = (m1 + m2 - 2 * k) / (x1 - x2) / (x1 - x2);
    b = ((m1 - m2) / (x1 - x2) - 3 * a * (x1 + x2)) / 2;
    c = m1 - 3 * a * x1 * x1 - 2 * b * x1;
    d = f1 - a * x1 * x1 * x1 - b * x1 * x1 - c * x1;
    v.push_back(a), v.push_back(b), v.push_back(c), v.push_back(d);
    return v;
}

// ppForm中可求M_i的情况求系数
vector<double> solve3(double x1, double x2, double f1, double f2, double M1, double M2)
{
    vector<double> v;
    double a, b, c, d;
    a = (M1 - M2) / (x1 - x2) / 6;
    b = (M1 - 6 * a * x1) / 2;
    c = (f1 - f2) / (x1 - x2) - a * (x1 * x1 + x1 * x2 + x2 * x2) - b * (x1 + x2);
    d = f1 - a * x1 * x1 * x1 - b * x1 * x1 - c * x1;
    v.push_back(a), v.push_back(b), v.push_back(c), v.push_back(d);
    return v;
}

// 把数组输入进向量
void pushback(vector<double> &x, double a[], int n)
{
    x.clear();
    for (int i = 0; i < n; i++)
        x.push_back(a[i]);
}
// 输出向量
void printvector(vector<double> v)
{
    for (int i = 0; i < v.size(); i++)
        cout << v[i] << " ";
    cout << endl;
}
// 求最大值
double maxvalue(vector<double> v)
{
    double x = v[0];
    for (int i = 1; i < v.size(); i++)
        if (v[i] > x)
            x = v[i];
    return x;
}

void Splines ::printspline()
{
    cout << x.size() << endl;
    printvector(x);
    for (int i = 0; i < coe.size(); i++)
        printvector(coe[i]);
}
void Splines ::fileprintspline(string s)
{
    fstream f;
    f.open(s, ios::out);
    f << x.size() << endl;
    for (int i = 0; i < x.size(); i++)
        f << x[i] << " ";
    f << endl;
    for (int i = 0; i < coe.size(); i++)
    {
        for (int j = 0; j < coe[i].size(); j++)
            f << coe[i][j] << " ";
        f << endl;
    }
    f.close();
}
double Splines ::value(double v)
{
    int n = coe[0].size();
    for (int i = 0; i < x.size(); i++)
    {
        if (v >= x[i] && v <= x[i + 1])
        {
            if (n == 2)
                return coe[i][0] * v + coe[i][1];
            else
                return coe[i][0] * v * v * v + coe[i][1] * v * v + coe[i][2] * v + coe[i][3];
        }
    }
    return 0;
}
#endif